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AN  ENERGY  DECAYING  SCHEME  FOR  NONLINEAR  DYNAMICS  OF  SHELLS 

CARLO  L.  BOTTASSO*,  OLIVIER  A.  BAUCHAUt,  AND  JOU-YOUNG  CHOI* 

Abstract.  A  novel  integration  scheme  for  nonlinear  dynamics  of  geometrically  exact  shells  is  developed 
based  on  the  inextensible  director  assumption.  The  new  algorithm  is  designed  so  as  to  imply  the  strict  decay 
of  the  system  total  mechanical  energy  at  each  time  step,  and  consequently  unconditional  stability  is  achieved 
in  the  nonlinear  regime.  Furthermore,  the  scheme  features  tunable  high  frequency  numerical  damping  and 
it  is  therefore  stiffly  accurate.  The  method  is  tested  for  a  finite  element  spatial  formulation  of  shells  based 
on  mixed  interpolations  of  strain  tensorial  components  and  on  a  two-parameter  representation  of  director 
rotations.  The  robustness  of  the  scheme  is  illustrated  with  the  help  of  numerical  examples. 

Key  words,  geometrically  exact  shell,  time  integration,  energy  decaying  scheme 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction.  The  formulation  of  integration  algorithms  for  nonlinear  dynamics  of  geometrically 
exact  shells  is  the  focus  of  this  work.  The  partial  differential  equations  governing  this  class  of  problems  are 
known  to  present  a  rich  mathematical  structure.  In  particular,  the  resulting  models  are  Hamiltonian  systems 
characterized  by  a  symplectic  nature  and  associated  with  conservation  laws  that  stem  from  symmetries  of 
the  Hamiltonian.  The  linear  and  angular  momentum  as  well  as  the  total  mechanical  energy  are  conserved 
for  free  motions  of  such  systems. 

The  understanding  of  the  geometric  characteristics  of  the  governing  equations  has  been  historically 
confined  to  the  fields  of  analytical  mechanics  and  pure  mathematics.  Surprisingly,  this  knowledge  has  been 
seldom  used  for  the  development  of  numerical  methods.  Indeed,  the  study  of  new  integration  algorithms  has 
been  traditionally  preoccupied  with  the  development  of  methods  applicable  to  vast  classes  of  problems,  for 
example  the  class  of  differential/algebraic  equations,  or  hyperbolic  conservation  laws.  Consequently,  classical 
methods  rarely  preserve  the  underlying  structure  of  the  problem  being  solved,  and  hence,  such  structure  is 
lost  in  the  numerical  solution. 

This  approach  also  limits  the  possible  theoretical  analyses  of  the  schemes,  which  are,  more  often  than 
not,  confined  to  linear  or  model  cases.  For  instance,  it  is  customary  to  characterize  integration  schemes 
for  structural  dynamics  by  studying  their  behavior  when  applied  to  a  linear  oscillator.  This  approach  is 
clearly  not  adequate  when  dealing  with  highly  nonlinear  problems  such  as  the  dynamics  of  geometrically 
exact  shells. 

A  new  approach  to  the  design  of  integration  algorithms  attempts  to  bridge  the  divide  between  theoretical 
and  numerical  mechanics.  Under  this  new  paradigm,  numerical  schemes  are  “backward-engineered”  to 
preserve  some  important  qualitative  features  of  the  governing  equations.  Fittingly,  this  approach  is  now 
called  geometric  integration  in  the  mathematical  community  [18].  Attempts  at  designing  “geometry-aware” 
algorithms  for  structural  dynamics  problems  can  be  traced  back  to  the  work  of  Simo  and  co-workers  who 
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analyzed  the  problems  of  rigid  body  dynamics  [26],  nonlinear  elasto-dynamics  [23],  geometrically  exact 
shells  [24],  and  geometrically  exact  beams  [25].  In  all  cases,  the  idea  was  to  design  algorithms  that  ensure 
the  discrete  preservation  of  the  total  mechanical  energy  and  of  the  linear  and  angular  momenta  of  the  system. 

When  integrating  linear  and  nonlinear  finite  element  models,  the  implications  of  the  discrete  equations 
stiffness  must  be  carefully  considered.  Indeed,  high  frequencies  are  an  artifact  of  the  spatial  discretization 
process  and  do  not  reflect  the  high  frequencies  of  the  original  infinite  dimensional  problem.  The  need  for 
high  frequency  numerical  dissipation  has  been  recognized  in  the  past  for  linear  problems  [20] .  When  dealing 
with  complex  nonlinear  systems,  numerical  dissipation  becomes  indispensable.  Indeed,  nonlinearities  provide 
a  mechanism  for  transferring  energy  from  low  to  high  frequency  modes.  Consequently,  numerical  solutions 
feature  violent  oscillations  of  a  purely  numerical  origin  that  will  eventually  play  havoc  with  the  convergence 
characteristics  of  the  nonlinear  equation  solver. 

Among  the  various  geometric  characteristics  of  shell  equations,  energy  preservation  appears  to  be  the 
most  important  for  the  development  of  robust  time  integration  schemes.  In  fact,  strict  energy  preservation 
at  the  discrete  level  leads  to  unconditional  stability  in  the  nonlinear  regime ,  whereas  the  classical  approach 
based  on  the  analysis  of  the  spectral  radius  leads  to  unconditional  stability  in  the  linear  regime  only.  An 
energy  preserving  (EP)  scheme  for  geometrically  exact  shell  is  developed  in  this  paper.  In  addition,  the 
scheme  also  preserves  both  linear  and  angular  momenta  of  the  system  at  the  discrete  level.  Unfortunately, 
preservation  of  energy  and  high  frequency  dissipation  cannot  coexist,  unless  energy  is  transferred  from  high 
to  low  frequency  modes,  a  transfer  that  has  no  physical  basis.  To  solve  this  problem,  a  family  of  energy 
decaying  (ED)  schemes  that  imply  a  controllable  energy  decay  within  each  time  step  is  proposed  in  this  work. 
In  geometric  terms,  this  means  that  the  evolution  of  the  system  is  not  confined  to  the  level  set  of  constant 
energy,  but  is  allowed  to  drift  away  from  it  in  a  monotonic  and  controllable  manner.  Since  the  energy  remains 
bounded  at  all  times,  the  scheme  is  unconditionally  stable  for  nonlinear  systems.  Furthermore,  it  can  be 
shown  that  the  energy  dissipation  mechanism  of  the  algorithm  is  the  result  of  the  removal  of  the  higher 
frequencies  from  the  computed  response. 

In  related  papers,  various  energy  preserving  and  decaying  geometric  integrators  were  developed  for  rigid 
bodies  and  geometrically  exact  beams  [10,  11,  14,  6,  7],  and  nonlinear  elastodynamics  [9].  The  concepts 
were  extended  to  multibody  systems  featuring  nonlinear  holonomic  constraints  [3,  15,  12,  16,  13];  non- 
holonomic  and  unilateral  constraints  were  treated  in  [5,  4].  The  integration  of  the  present  shell  model  in  a 
general  finite  element  based  multibody  framework  is  discussed  in  [8].  The  proposed  scheme  is  independent 
of  the  choice  of  spatial  discretization  applied  to  the  governing  partial  differential  equations.  In  the  present 
implementation,  the  finite  element  method  is  used,  and  the  mixed  interpolation  of  tensorial  components  [1, 
2,  17]  is  implemented  to  avoid  the  shear  locking  problem.  The  orientation  of  unit  shell  directors  is  described 
by  a  special  family  of  two-parameter  rotations. 

The  paper  is  laid  out  as  follows.  The  classical  equations  of  motion  for  geometrically  exact  shells  based 
on  inextensible  unit  directors  are  presented  in  section  2.  Next,  an  EP  scheme  is  developed  in  section  3. 
Section  4  then  presents  an  ED  algorithm  with  tunable  high  frequency  dissipation  that  is  constructed  from 
the  EP  scheme.  Finally,  numerical  examples  are  presented  in  section  5  to  demonstrate  the  efficiency  and 
robustness  of  the  proposed  scheme.  A  discussion  section  concludes  the  paper. 

2.  Formulation  of  the  Equations  of  Motion. 

2.1.  Kinematics  of  the  Shell  Problem.  Consider  a  shell  of  thickness  h  and  reference  surface  area  Q, 
as  depicted  in  fig.  2.1.  An  inertial  frame  of  reference  S  consisting  of  three  mutually  orthogonal  unit  vectors 
ii,  i2,  i3  is  used.  Let  be  the  position  vector  of  an  arbitrary  point  on  the  reference  surface  of  the  shell,  and 
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let  £  be  the  material  coordinate  along  n,  the  normal  to  the  reference  surface.  The  position  vector  r  of  an 
arbitrary  point  on  the  shell  in  its  reference  configuration  is  then 

r(C,£2. 0  =£o(C>£2)  +  C«(C,£2).  C2-1) 


where  £*and  £2  are  the  material  coordinates  used  to  represent  the  shell  reference  surface.  The  coordinates 
}  £2  and  £  form  a  set  of  curvilinear  coordinates  that  are  a  natural  choice  to  represent  the  shell  geometry. 
The  coordinates  £*and  £2  are  assumed  to  be  lines  of  curvatures  of  the  shell  reference  surface.  The  base 
vectors  are  then 


9  = 


r  i 

dr 

ll 

E.2.  af 

(1 —  ^r)  — i ^  (1  )  fia*  4 


(2.2) 


where  R\  and  R2  are  the  principal  radii  of  curvature,  aa  =  and  the  notation  (*)>Q  is  used  to  denote  a 
derivative  with  respect  to  .  It  is  convenient  to  introduce  a  set  of  three  mutually  orthogonal  unit  vectors 
at  the  shell  reference  surface  ( i.e .  at  C  —  0) 


Ci  = 


4i  . 
v^n  ’ 


£3  =  n, 


(2.3) 


where  aaa  =  a a  *  aa. 

Two  fundamental  assumptions  will  be  made  concerning  the  deformation  of  the  shell,  i.e.  the  material 
line  initially  normal  to  the  reference  surface  of  the  shell  remains  a  straight  line  and  suffers  no  extension. 
This  is  the  classical  inextensible  director  model.  With  these  assumptions,  the  position  vector  of  a  material 
point  of  the  shell  writes 


m1 ,  2 !  o = ro  1 ,  e ) + uk1 


(2.4) 
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where  u(^ .  £2)  is  the  reference  surface  displacement  vector.  In  the  deformed  configuration,  the  base  vectors 
at  the  shell  reference  surface  are 


G  =  [£,,  G2,  £3]  = 


d&\ 

0CJ ' 


(2.5) 


Introducing  the  position  vector,  eq.  (2.4),  then  yields 


fil  &  r  ' 

V5IT  y/a£’  -3 


=  E  +  CH, 


(2.6) 


where 


E  —  [£i ?  £2?  £3]  — 


U,2 


-1  +  v^t’  S2  +  vW  ’  “  Lvsir  %/®22 1 

Note  that  £3(^5  £2)  is  a  unit  vector,  whereas  Ej  and  £2  are  not  unit  vectors,  nor  are  they  orthogonal  to 
IS3,  as  axial  and  transverse  shearing  strains  develop  during  deformation. 


H  = 


£3,1  £3,2 


,0 


(2.7) 


2.2.  Equations  of  Motion.  The  Green-Lagrange  strain  tensor  e  is  defined  as 

e=  ^(GTG-gTg).  (2.8) 

The  strain  tensor  e  is  defined  in  the  curvilinear  coordinate  system  defined  by  coordinates  £1,£2  and  (. 
However,  it  is  more  convenient  to  work  with  the  strain  tensor  e  defined  in  the  locally  rectangular  system 
defined  by  triad  ex,  e2,  £3,  see  eqs.  (2.3).  For  shallow  shells  ( i.e .  C/Ri  C  1  and  C/R2  1)  undergoing  large 
displacements  and  rotations  but  small  strains  (all  strain  components  are  assumed  to  be  small  compared  to 
unity),  the  strain-displacement  relationships  can  be  written  as 


where 


e  =  J  [EtE  - 1  +  C  (. EtH  +  HtE  +  «)]  , 


l/i?i  0  0 

0  1/R2  0 

0  0  0 


(2.9) 


(2.10) 


It  is  clear  that  the  strains  can  be  expressed  in  terms  of  five  parameters:  the  three  components  of  the 
displacement  field  u  (through  Et  and  E2)  and  the  two  parameters  defining  the  orientation  of  the  unit 
director  e3.  Virtual  changes  in  the  strain  energy  of  the  structure  are  given  by 

6V=  f  f  SV  d(dn=  f  [  6e-rd(dn,  (2.11) 

Jn  Jh  Jn  Jh 

where  8V  is  the  virtual  strain  energy  density,  and  r  the  second  Piola-Kirchhoff  stress  tensor.  Introducing 
the  strains,  eq.  (2.9),  and  taking  into  account  the  symmetry  of  the  stress  tensor  then  yields 


8V  =  SE  •  (E  +  (H)t  +  8H  ■  QEr. 


(2.12) 


The  existence  of  a  strain  energy  density  function  V  is  postulated  here,  hence  the  constitutive  laws  are  of  the 
form  r  —  dV  jde . 

The  velocity  vector  of  material  point  P  of  the  shell  is  obtained  by  differentiating  the  position  vector, 
eq.  (2.4),  with  respect  to  time,  to  find  v  =  u  +  (£3-  The  kinetic  energy  of  the  system  is  now 

K  =  [  [  K  d(dft  =  i  f  f  pv-v  d(dH,  (2.13) 

Jn  Jh  2  Jq  Jh 
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where  K  is  the  kinetic  energy  density.  Introducing  the  velocity  vector  then  yields 

^=^p(m  +  c4)-(m  +  C4)- 

Hamilton’s  Principle  can  now  be  expressed  as 

ftf  [  f  (SK  —  6V)  dCdftd t 

Jti  J n  Jh. 

=  f'  f  Jh  [p(S M  +  C$Es)  •  (u  +  (Ms)  +  SE  ■  (E  +  (H)r  +  SH  ■  ( Et]  dCdOdi 


(2.14) 


=  0.  (2.15) 


Integrating  through  the  thickness  of  the  shell,  we  get 

J  J  *  [A  —  {£Li,i  +  jy2|g)]  +  SE_  ■  ^  —  (Ml, I  +  M2, 2)  “b  £3]  |  dfldt  =  0.  (2.16) 

In  this  expression,  h  —  mu  4-  s^E^,  and  g  =  s*u  -f  i*k 3  are  the  linear  and  angular  momentum  vectors  of 
the  shell,  respectively;  the  mass  coefficients  are  defined  as  m  —  Jhp  d£,  s*  =  fh  p(  d£,  I*  =  fhP(2 
The  in-plane  forces  are  N_a  =  (EjV*  +  HM*a)/^/a„r*<  the  out-of-plane  forces  jVg  —  ETQ,  and  the  bending 
moments  Ma  =  (EM^/^/aM*  The  convected  forces  are  N*  =  [iV*t  JVg^IVg]  =  fhr  d(,  and  the  convected 
bending  moments  M*  =  [Mi 5 M2? M3]  =  fhTC 

The  equations  of  motion  of  shells  could  be  derived  from  this  principle  by  expressing  the  variations  SE 3 
in  terms  of  two  components  of  virtual  rotation. 

3.  Energy  Preserving  Scheme.  Discrete  equations  of  motion  that  imply  discrete  conservation  laws 
for  the  total  mechanical  energy,  linear  momentum  and  angular  momentum  of  the  system  will  now  be  de¬ 
veloped.  Times  tj  and  tf  denote  the  initial  and  final  times  for  a  time  step,  respectively,  and  the  subscripts 
(•)*  and  (• )f  indicate  quantities  at  U  and  tf,  respectively.  Furthermore,  the  subscript  (-)m  is  used  to  denote 
mid-point  average  quantities  defined  as 

(0m  =  |  [(■)/  +  (•)<]•  (3-1) 

The  following  matrix  identity  will  be  used  extensively 

AjBf-AjBi  =  (Af-Ai)TBm+Al(Bf-Bi).  (3.2) 

Hamilton’s  Principle,  eq.  (2.15),  is  now  approximated  in  time  in  the  the  following  manner 


/o  Jh 


P  [(«/  “  Hj)  +  C(^3/  —  — 3i)] 


Uf -Mi 


+c 


£3/  ~  Ejdi 


+(Ef  -  Ei)  ■  (Em  +  (Hm)Ta  +  (Hf  -  H^)  ■  CEmTa  j  dCdft  =  0.  (3.3) 

The  change  in  strain  components  from  to  tf  is  evaluated  with  the  help  of  identity  (3.2)  to  find 


ef  -  ^  i  \(Ef  -  Ei)T(Em  +  (Hm)  +  (Em  +  C,Hm)T(Ef  -  E{) 
+(Hf  -  Hi)T(Em  +  (El(Hf  -  Hi)] . 


(3.4) 
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Over  one  time  step,  the  strain  components  can  be  approximated  as  e(rj)  =  em  4-  rj(ef  —  ei)j 2,  where  77  = 
2 (t  -  tm)l At  is  the  non-dimensional  time.  If  the  strain  energy  density  function  V  is  viewed  as  a  function  of 
the  scalar  variable  77,  the  mean  value  theorem  then  implies  the  existence  of  a  fj  6  [—1, 1]  such  that 


de 

dr] 


2  =  Vi +ra  •  (ef  -e<). 


(3.5) 


This  relationship  defines  the  average  second  Piola-Kirchhoff  stress  tensor,  ra  =  dV/de\^.  Combining  this 
result  with  eq.  (3.4)  then  leads  to 


(Ef  -  Ei)  ■  (Em  +  C Hm)Ta  +  ( Hf  -  Hi)  ■  C EmTa  =  (e,  -  a)  ■  ra  =  V,  -  %  (3.6) 


where  the  symmetry  of  the  stress  tensor  was  taken  into  account.  For  linear  constitutive  laws  of  the  form 
r  —  C*  e.  where  C*  is  the  stiffness  matrix,  the  average  stress  tensor  simply  becomes  ra  =  C*  em. 

The  following  configuration  updates  are  now  defined 


Uf  —  u2 

At 


£3/  ~3i  _  p 

At  ~ 3m' 


(3.7) 


Introducing  eqs.  (3.6)  and  (3.7)  into  the  approximate  expression  for  Hamilton’s  Principle,  eq.  (3.3),  then 
leads  to 


I  J  { f  (Urn  +  CE3m)  •  [(uf  -  Ui)  +  C(E3f  -  4i)]  +  (Vf  -  Vi) }  d<dO 

=  L  L  {  2  ^  +  C^f) '  +  C^f)  ~2{ki  +  C^i] '  +  C-3i)  +  ~  fi)}  dcdn 

=  [  f  [( Kf  -  Ki)  +  (Vf  -  Vi)]  d(dQ  =  0.  (3.8) 
Jn  Jh 


This  result  clearly  implies  the  conservation  of  the  total  mechanical  energy  of  the  system  within  a  step. 
In  summary,  the  approximate  form  of  Hamilton’s  Principle  given  by  eq.  (3.3)  leads  to  a  discrete  energy 
conservation  statement,  eq.  (3.8),  when  the  configuration  updates  are  chosen  according  to  eqs.  (3.7),  and  the 
average  stress  according  to  eq.  (3.5). 

Integrating  through  the  thickness  of  the  shell  leads  to 


^— ^ -(&«.!+ JSam.a) 


At 


+  ( E3f  -  E^i)  ■ 


if-ii 


At 


(M-lm,l  +  M.2m,2)  +  i^3mj  j 


dfl  =  0. 


(3.9) 


In  this  expression,  the  in-plane  forces  are  Nam  =  (EmN_*aa  -f  HmM_*aa)/ ^/aaa,  the  out-of-plane  forces  jV3m  — 
EmKL,  and  the  bending  moments  Mam  =  (£mM)m)/v/<w  The  discrete  governing  equations  of  motion 
for  shells  are  then 

^  -«,„..+*£*»)  =&,;  M 

Qm  ~  +  =  C~  <3  U> 

where  p  are  the  externally  applied  loads,  and  q*  the  externally  applied  moments  measured  in  the  local 
system.  The  finite  change  in  director  orientation  —  Esi  was  expressed  in  terms  of  the  two  parameter 
incremental  rotation  vector,  see  B. 
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Invariance  of  the  system  Hamiltonian  under  spatial  translations  and  rotations  implies  the  conservation 
of  the  linear  and  angular  momenta.  Although  discrete  preservation  of  momenta  is  less  crucial  than  discrete 
preservation  of  energy,  it  is  interesting  to  note  that  eqs.  (3.10)  and  (3.11)  also  imply  the  discrete  preservation 
of  this  invariant.  At  first,  eqs.  (3.10)  are  projected  onto  the  test  functions  5r  (r^  -f  um)  and  eqs.  (3.11)  onto 
the  test  functions  4  £3 m,  where  7r  is  an  arbitrary  vector.  Next,  integration  over  the  shell  reference  surface 
yields 


Ja  jiffeo  +  Um)  • 


h}  -hj 


At 


~  (— lm,l  +  N_2ma) 


'If -Si 

At 


1  -  (Kim, l  +  M2m,2)  +  Ksm  }  dfl  =  0.  (3. 


12) 


Straightforward  algebraic  manipulations  then  lead  to 


^  •  J  [go  +  uf)  hf  -  (eq  +  Hi)  hi  +  hm(uf  -  Mi) 

+  E.f9_f  ~  E-igi  +  iJEf  -  i?j)]  dQ  =  0,  (3.13) 

where  the  following  result  was  used 

E-lmH-lm  —2m— 2m  T  +  Ez,2mM-.2m  +  Esmitsm  “  (3.14) 


Inserting  the  configuration  updates,  eqs.  (3.7),  into  eq.  (3.13)  then  yields 

ft' f[fc 0  +  2/)  hf-(r o  +  Mi)  hi  +  E3fgf  -  +  (hmum  +  §m  £3™)]  dtt  =  0.  (3.15) 

It  is  easily  verified  that  hmUm  +  gmEsm  -  0.  Hence,  since  tt  is  arbitrary,  eq.  (3.15)  implies  the  discrete 
conservation  of  the  total  angular  momentum,  Jfi[(r0  +5)  k  +  £L  £]  dD.  Finally,  projecting  eqs.  (3.10)  onto 
the  test  functions  7r  and  eqs.  (3.11)  onto  the  null  test  functions  gives  the  discrete  conservation  of  the  total 
linear  momentum  fQ  h  dH. 

It  is  important  to  note  that  any  spatial  discretization  of  the  discrete  equations  of  motion  will  inherit 
the  discrete  energy  and  momentum  conservation  statements  just  proved,  when  the  configuration  updates  are 
chosen  according  to  eqs.  (3.7),  and  the  average  stress  according  to  eq.  (3.5). 


4.  Energy  Decaying  Scheme.  As  discussed  in  the  introduction,  energy  preservation,  per  se ,  is  not 
sufficient  to  yield  robust  time  integration  schemes.  High  frequency  numerical  dissipation  must  be  added  as 
an  inherent  feature  of  the  scheme.  Such  a  scheme  will  now  be  constructed  for  the  shell  equations  of  motion 
using  the  EP  scheme  as  a  basic  building  block. 

First,  an  additional  state  is  introduced  at  time  tj  =  \ime-+o(ti  4*  e),  and  the  subscript  (-)j  is  used  to 
denote  quantities  at  this  time.  The  following  averages  are  now  defined 

(•)*  =  5  [(■)/  + (•)>];  (•)*  =  5  [(•)*  + (•)«]•  (4-1) 

The  ED  scheme  proceeds  from  the  initial  to  the  final  time  by  means  of  two  coupled  steps:  one  step  from  ti 
to  tf,  the  other  from  U  to  tj.  The  time-discrete  equations  of  dynamic  equilibrium  are 


Q 


^  +^29,2)  ~Pg’ 

T  T  ~  Q +  ^9,2  -  K3g)  =  g*; 


(4.2) 
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(4.3) 


\t  +  5  +  ^23,2)  -  (^ip,i  +  ^,2)]  -  £h; 


Qh  J \+  +  o  IPs  (Mi9>1  +i^29,2  ~  H-3g)  ~  Qh  (Mlp.l  +  M2p,2  i^p)]  ~  2ft  ‘ 


The  configuration  update  relationships  are  given  as 

Uf=Ui  +  At  ( uf  +  Uj)/ 2,  Uj  =  u{  -  At  [uf  -  Uj  -  a(Uj  -  «j)]  /6; 

£3/  =  £34  +  At  (4/  +  iw)/2,  =Esi- At  [4,  -  4,  -  «(4i  -  £w)]  /6. 


(4.4) 


where  a  is  a  tuning  parameter  that  controls  the  amount  of  numerical  dissipation  provided  by  the  scheme, 
while  the  forces  7VQp  and  moments  M_ap  are  given  by 


Kap  —  N_ah  +  a  (Naj  -  Nai)/ 2;  M,p  =  +  a  (M*  "  Mai)/2- 


(4.5) 


Using  developments  similar  to  those  exposed  for  the  EP  scheme,  it  can  be  easily  shown  that  the  proposed 
discrete  equations  imply 


(Kf  +  Vf)  -  (Ki  +  Vi)  +  ac2=  0. 


(4.6) 


c2  is  a  positive  quantity  given  by 

c2  =  Jnl  [m II it II  •  11  i II  +2s*  II i II  •  II 4 II  +r  || 4 II  •  II 4  ll] dn 

+  /  5  ||  fi  ||  c*  |U  II  dn  >  0,  (4-7) 

where  ||  •  ||=  (-)j  —  (•)*  is  the  jump  between  ti  and  tj.  This  result  implies  the  decay  of  the  total  mechanical 
energy  over  one  step  of  the  algorithm,  (Kf  +  Vf)  <  (Ki  +  Vi).  The  parameter  a  clearly  controls  the  amount 
of  energy  that  is  dissipated  within  the  step.  Two  such  parameters  could  be  used,  controlling  the  amount  of 
dissipated  kinetic  and  strain  energies,  respectively,  but  this  level  of  complexity  does  not  seem  to  be  necessary. 
The  property  of  preservation  of  momentum  observed  in  the  EP  case  is  lost  in  the  ED  algorithm. 

If  the  above  ED  scheme  is  applied  to  a  single  degree  of  freedom  linear  oscillator,  the  asymptotic  value 
of  the  spectral  radius  of  the  amplification  matrix,  poo,  is  found  to  be  poo  =  (1  —  Ck)/(1  +  a).  For  a  =  1, 
Poo  —  0,  and  asymptotic  annihilation  is  achieved.  If  a  =  0,  poo  =  1?  and  in  view  of  eq.  (4.6),  energy  is 
exactly  preserved.  Hence,  the  ED  scheme  is  in  fact  a  family  of  schemes  with  a  single  tuning  parameter, 
a,  that  controls  the  amount  of  high  frequency  numerical  dissipation;  both  asymptotic  annihilation  or  exact 
energy  preservation  can  be  achieved  with  the  same  scheme  by  using  a  =  1  or  0,  respectively. 

5.  Numerical  Examples.  All  the  examples  described  in  this  section  will  be  treated  with  the  proposed 
ED  family  of  schemes  corresponding  to  values  of  the  tuning  parameter  a  €  [0,1].  Although  any  value  of  a 
within  this  range  can  be  used,  the  examples  described  here  will  contrast  the  two  extreme  choices.  For  a  =  1 
(p^  =  0),  asymptotic  annihilation  is  obtained,  and  this  will  be  called  the  ED  scheme.  On  the  other  hand, 
for  a  —  0  (poo  =  1),  exact  energy  preservation  is  achieved,  and  this  will  be  called  the  EP  scheme. 

5.1.  Clamped  Half-Cylinder  under  Point  Load.  Consider  the  half-cylinder  of  radius  R  =  1.2  m 
and  width  b  =  2  m  depicted  in  fig.  5.1.  The  shell  has  a  thickness  t  =  6  mm,  is  build-in  along  edge  BC 
and  free  along  the  other.  The  structure  is  made  of  aluminum;  Young’s  modulus  E  =  73  GPa,  Poisson’s 
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Fig.  5.2.  Configuration  of  the  system  at  various  instants  in  time. 


ratio  v  =  0.30  and  density  p  =  2700  kg/m3.  At  point  D,  the  shell  is  subjected  to  a  concentrated  load 
P  —  —  Po(t)(i[  +i2  +  is).  The  magnitude  of  the  load  is 


Po(t)  = 


P  (l-cos27rt/T)/2  t  <T, 
0  t>T, 


(5.1) 


where  P  =  0.1  kN  and  T  =  2.0  s.  The  shell  was  modeled  by  a  regular  8x4  mesh  of  quadratic  elements.  All 
simulations  were  run  with  a  time  step  At  =  5.0  10-03  s,  for  a  total  time  of  6  s. 

Under  the  effect  of  the  applied  loads,  the  shell  bends  predominantly  in  the  vertical  direction,  its  direction 
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Fig.  5.3.  Displacement  components  at  point  D.  Ui :  A;  U2 :  o;  U$:  □.  ED  scheme:  solid  line;  EP  scheme:  dashed  line. 


Fig.  5.4.  Time  history  of  bending  moments  at  point  M.  M\\:  A;  M12:  o;  M22;  □  .  ED  scheme:  solid  line ;  EP  scheme: 
dashed  line. 

of  least  bending  stiffness,  as  illustrated  in  fig.  5.2  that  shows  the  configuration  of  the  system  at  various  instants 
in  time.  The  three  components  of  displacements  at  point  D  are  shown  in  fig.  5.3;  vertical  displacements 
of  up  to  0.6  m  are  observed.  The  bending  and  twisting  moments,  Mu,  M2 2,  and  M12,  respectively,  at 
point  M  are  shown  in  fig.  5.4.  Note  the  significant  transverse  and  twisting  moments  associated  with  the 
three-dimensional  motion  of  the  shell.  The  components  of  in-plane  and  transverse  shearing  forces  are  shown 
in  fig.  5.5  and  5.6,  respectively. 

Next,  the  same  problem  was  simulated  with  p =  1,  Le.  with  no  high  frequency  dissipation.  The 
corresponding  results  are  shown  in  figs.  5.3  to  5.6.  Displacement  and  moment  results  are  found  to  be  in 
excellent  agreement.  At  the  scale  of  the  figures,  they  are,  in  fact,  indistinguishable.  For  the  period  of  time 
2  <  t  <  6  s,  the  system  is  not  subjected  to  any  loading,  and  the  total  mechanical  energy  of  the  system  should 
remain  constant.  For  the  EP  scheme,  the  energy  is  indeed  preserved,  as  expected;  for  the  ED  scheme,  0.3% 
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Fig.  5.5.  Time  history  of  in-plane  forces  at  point  M.  Fn:  top  figure;  F12 :  middle  figure;  F22 ■  bottom  figure.  ED  scheme: 
solid  line;  EP  scheme:  dashed  line. 


Fig.  5.6.  Time  history  of  transverse  shear  forces  at  point  M.  F13;  top  figure;  ^23*  bottom  figure.  ED  scheme:  solid  line ; 
EP  scheme:  dashed  line. 


of  the  energy  is  numerically  dissipated  in  this  period.  It  could  be  concluded  that  the  EP  and  ED  solutions 
are  nearly  identical,  and  that  numerical  dissipation  is  not  necessary.  However,  the  ED  and  EP  scheme 
predictions  for  the  in-plane  and  transverse  shearing  forces,  shown  in  fig.  5.5  and  5.6,  are  markedly  different. 
EP  predictions  for  force  components  Fn,  F12,  and  F13  show  high  frequency  oscillations  that  are  absent  in 
the  corresponding  ED  predictions.  A  simulation  using  the  ED  scheme  with  a  time  step  At  =  1.0  10  03  s 
showed  that  the  ED  predictions  are  converged.  A  simulation  using  the  EP  scheme  and  the  same  smaller  time 
step  yielded  results  with  increased  high  frequency  oscillations  for  the  force  predictions.  It  should  be  noted 
that  the  dynamic  response  of  this  simple  system  is  very  smooth;  yet  even  here,  high  frequency  numerical 
dissipation  appears  to  be  necessary  to  obtain  a  smooth,  converged  solution. 
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b=0.05m 


Fig.  5.7.  Configuration  of  the  plate  with  concentrated  masses. 


Fig.  5.8.  Time  history  of  the  total  mechanical  energy  (top  figure)  and  relative  energy  loss  (bottom  figure).  Plate  model 
(EDS):  solid  line;  Plate  model  (EPS):  dashed  line;  Beam  model  (EDS):  dashed-dot  line . 


5.2.  Dynamic  Response  of  a  Plate  with  Edge  Beams.  Consider  the  rectangular  plate  of  length 
L  =  2  m,  width  b  =  0.05  m  and  thickness  t  =  2  mm  as  depicted  in  fig.  5.7.  Two  circular  beams  of  radius 
r  =  2  mm  are  attached  at  the  plate  edges.  A  third  circular  beam  is  located  at  the  center  of  the  plate.  All 
components  are  made  of  aluminum;  Young’s  modulus  E  =  73  GPa,  density  p  =  2700  kg/m3.  The  total  mass 
of  the  edge  beams  is  10  kg  each,  and  that  of  the  central  beam  is  1  kg.  The  plate  is  subjected  to  uniformly 
distributed  loads  Fm  and  Ft  along  FE  and  CB,  respectively.  The  components  of  these  loads  along  the  ix 
and  i3  axes  are  F\m  =  40  N/m  and  F3m  —  80  N/m,  respectively,  and  Flt  =  -20  N/m  and  Fu  =  -60  N/m, 
respectively.  The  common  time  history  of  each  loading  component  is 


Fi(t)  = 


Fi  (1  -cos27rf/T)/2  t  <T, 
0  t>T, 


(5.2) 
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Fig.  5.9.  Trajectory  of  the  plate  mid-span  point.  Plate  model  (EDS):  solid  line;  Plate  model  (EPS):  dashed  line;  Beam 
model  (EDS):  dashed-dot  line. 


Fig.  5.10.  Time  history  of  the  quarter-span  axial  force .  Plate  model  (EDS):  solid  line;  Plate  model  (EPS):  dashed  line; 
Beam  model  (EDS):  dashed-dot  line. 


where  T  =  3.0  s.  The  plate  is  discretized  with  10  quadratic  plate  elements  along  its  length.  A  constant  time 
step  At  —  6  10~03  s  was  used  for  all  simulations. 

At  time  t  >  3  s  the  applied  load  vanishes,  and  the  system  total  mechanical  energy  should  remain 
constant.  The  evolution  of  the  energy  of  the  system  for  three  different  cases  is  shown  in  fig.  5.8:  the  plate 
model  using  the  ED  scheme,  the  same  plate  model  using  the  EP  scheme,  and  a  simplified  model  of  the  system 
using  beam  elements  and  the  ED  scheme.  For  times  t  >  3  s,  the  total  energy  remains  a  constant  for  the  EP 
scheme  and  nearly  constant  for  the  ED  schemes.  The  relative  energy  loss  is  also  presented  in  the  figure:  0.6% 
of  the  energy  was  dissipated  by  the  ED  scheme  in  the  period  t  €  [3, 20]  s.  This  figure  clearly  demonstrates 
the  non-increasing  property  of  the  energy  evolution  for  the  proposed  ED  schemes,  as  opposed  the  constant 
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Fig.  5.11.  Time  history  of  the  quarter-span  transverse  shear  force .  Plate  model  (EDS):  solid  line;  Plate  model  (EPS): 
dashed  line;  Beam  model  (EDS):  dashed- dot  line . 


energy  predicted  in  EP  simulations.  The  trajectory  of  the  plate  mid-span  point  is  shown  in  fig.  5.9:  good 
correlation  is  observed  between  the  predictions  of  the  three  models.  The  beam  model  is  slightly  off  due  to  the 
inherent  simplifying  assumptions.  The  behavior  of  the  quarter- span  axial  force  and  transverse  shear  force 
are  shown  in  fig.  5.10  and  5.11,  respectively.  The  poor  predictions  of  the  EP  schemes  are  obvious  in  these 
two  plots.  The  history  of  axial  force  presents  violent  oscillations  with  amplitudes  an  order  of  magnitude 
larger  than  those  observed  for  the  ED  scheme.  The  history  of  the  transverse  shear  force  predicted  by  the 
EP  scheme  quickly  diverges  from  the  ED  predictions  for  both  beam  and  plate  models.  To  ascertain  the 
accuracy  of  the  ED  predictions,  a  convergence  study  was  performed.  Nearly  identical  results  were  found 
with  smaller  time  step  sizes  At  =  3.0  and  1.0  10“03  s,  or  when  using  the  time  adaptivity  procedure.  On 
the  other  hand,  oscillations  of  increasing  amplitude  were  found  as  the  time  step  size  is  reduced  in  the  EP 
scheme.  Furthermore,  the  time  adaptivity  procedure  failed  to  yield  any  results  because  the  time  step  size 
was  driven  to  unreasonably  small  values,  At  =  lO-07  s,  as  the  procedure  tries  to  cope  with  increasingly 
violent  oscillations. 
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Fig.  5.12.  Configuration  of  the  cruciform  problem. 


Fig.  5.13.  Configuration  of  the  cruciform  at  times  t  =  0.062  and  0.093  s. 


5.3.  Dynamic  Response  of  a  Cruciform.  Consider  a  cruciform  consisting  of  four  thin  panels  ( Panels 
A,  B,  C,  and  D )  connected  to  a  central  beam,  as  depicted  in  fig.  5.12.  Each  panel  is  of  thickness  t  =  4  mm, 
length  L  =  1.2  m,  and  width  b  =  0.1  m.  The  central  beam  has  a  square  cross-section  of  width  a  =  8  mm. 
A  mass  M  =  12  kg  is  attached  at  the  tip  of  the  central  beam  at  point  T.  Panels  and  beam  are  simply 
supported  at  the  root  of  the  cruciform.  A  concentrated  load  P(t)  is  applied  at  point  T .  The  load  acts  in 
the  plane  defined  by  axes  i2  and  £3  and  makes  a  30  degree  angle  with  axis  £2.  All  components  are  made  of 
aluminum  with  properties  given  in  the  previous  example.  The  time  history  of  the  applied  load  is 


P(t)  = 


P0  (1  —  cos27rt/T)/2  t  <  T, 
0  t>T, 


(5.3) 


where  P0  —  1-2  kN  and  T  =  0.1  s. 

As  the  applied  load  increases,  in-plane  stresses  in  the  panels  rapidly  increase  and  buckling  takes  place 
in  those  panels  subjected  to  compression,  as  can  be  observed  in  fig.  5.13  that  depicts  the  configuration  of 
the  cruciform  at  two  instants  in  time.  The  trajectory  of  point  T  projected  onto  plane  i2,  i3  is  shown  in 
fig.  5.14.  For  reference,  the  corresponding  trajectory  of  a  beam  with  cross-sectional  properties  equivalent  to 
those  of  the  cruciform  is  also  presented.  Of  course,  the  equivalent  beam  model  is  much  stiffer  since  it  does 
not  allow  buckling  to  take  place.  Furthermore,  the  motion  remains  confined  to  the  plane  defined  by  axis 
and  the  line  of  action  of  the  applied  load.  When  each  panel  is  modeled  individually,  the  stiffness  of  system 
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Fig.  5.14.  Trajectory  of  point  T  of  the  cruciform  projected  in  the  plane  is-  Solid  line:  shell  model;  dashed  line:  beam 
model. 


Fig.  5.15.  Total  mechanical  energy  of  the  system.  Solid  line:  shell  model;  dashed  line:  beam  model. 

varies  both  spatially  and  temporally,  giving  rise  to  the  more  complex  motion  shown  in  fig.  5.14.  The  total 
mechanical  energy  of  the  system  is  shown  in  fig.  5.15.  Prom  time  t  =  0.1  to  0.2  s,  the  system  is  free  and 
its  total  mechanical  energy  should  remain  constant.  Due  the  dissipative  nature  of  the  integration  scheme,  a 
small  amount  of  energy  is  dissipated  over  that  period  of  time:  2.7%  of  the  energy  was  dissipated  over  the 
2435  time  step  period. 

The  root  shear  force  and  quarter-span  bending  moment  in  Panels  A  and  C  are  shown  in  fig.  5.16 
and  5.17,  respectively.  Each  panel  undergoes  alternating  phases  of  tensile  and  compressive  loading.  During 
the  compressive  phases,  buckling  takes  place,  and  large  shear  forces  and  bending  moments  are  observed  in 
contrast  with  the  tensile  phases  during  which  these  quantities  remain  much  smaller. 
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Fig.  5.19.  System  configurations  at  various  time  instants  during  the  simulation. 


5.4.  Snap-Through  of  a  Cylindrical  Shell.  The  snap-through  behavior  of  a  cylindrical  shell  under 
a  concentrated  load  was  investigated  in  ref  [22].  The  shell  consists  of  a  60  degree  sector  of  a  cylinder  of 
height  h  =  5  m,  radius  R  =  5  m  and  thickness  t  =  0.1  m,  as  shown  in  fig.  5.18.  Material  properties  are: 
Young’s  modulus  E  =  210  GPa,  Poisson  ratio  v  =  0.25  and  density  p  =  104  kg/m3.  The  two  straight  edges 
of  the  shell  are  simply  supported,  while  the  two  curved  edges  are  free. 

A  concentrated  force  F  is  applied  at  the  shell’s  apex.  This  force  linearly  increases  from  0  to  5  107  N  in 
0.2  s,  then  is  held  constant  at  that  value.  The  simulation  ends  at  time  t  =  0.3  s.  Due  to  the  symmetry  of 
the  problem,  a  quarter  shell  only  is  modeled;  a  regular  4x4  mesh  of  quadratic  elements  was  used.  The  time 
step  size  was  selected  as  At  —  10“ 3  s. 

As  the  load  increases,  the  shell  apex  displacement  increases,  then  suddenly,  snap-through  takes  place 
and  curvature  reverses.  Curvature  reversal  initiates  in  the  region  of  the  applied  load,  then  quickly  propagates 
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Fig.  5.20.  Time  history  of  the  plate  center  vertical  displacement.  ED  scheme:  solid  line;  EP  scheme:  dashed  line. 


FlG.  5.21.  Time  history  of  the  plate  center  forces  for  the  ED  and  EP  schemes.  Force  Fx:  solid  line ;  F3:  dashed  line.  For 
clarity  the  EP  results  are  shifted  downwards  3  108  N. 


throughout  the  entire  structure,  which  undergoes  subsequent  violent  oscillations.  Snapshots  of  the  system 
at  various  instants  in  time  are  given  in  fig.  5.19.  The  vertical  displacements  of  the  point  of  application  of  the 
load  computed  with  the  ED  scheme  is  shown  in  fig.  5.20.  Note  the  gradual  increase  of  the  shell  deflection, 
until  collapse  at  buckling  and  the  resulting  vibratory  response  in  the  inverted  configuration. 

Ref.  [22]  presents  simulations  of  this  problem  using  various  schemes:  the  generalized-a  [19]  and  the 
CEMA  [22]  schemes.  The  former  scheme  features  high  frequency  numerical  dissipation  and  linear  stability 
properties,  while  the  latter  adds  to  the  generalized-cn  method  a  constraint  on  the  total  mechanical  energy  of 
the  system.  CEMA  is  therefore  both  energy  preserving  and  high  frequency  dissipative.  The  results  presented 
in  fig.  5.20  are  in  close  agreement  with  those  obtained  with  the  generalized-a  scheme,  but  quite  different 


19 


Fig.  5.22.  Time  history  of  the  plate  center  vertical  velocity.  ED  scheme:  solid  line;  EP  scheme :  dashed  line,  shifted 
downwards  3  102  m/s  for  clarity. 


from  those  predicted  by  CEMA  in  the  post  buckling  regime.  It  is  important  to  realize  that  the  higher 
modes  are  only  an  artifact  of  the  discretization  process,  and  should  therefore  be  removed  from  the  computed 
response.  A  standard  scheme  like  the  generalized-o  method  accomplishes  this  goal  through  the  characteristic 
low-pass  shape  of  its  spectral  radius;  however,  there  is  no  guarantee  that  energy  will  not  be  allowed  to  grow 
within  one  step  for  nonlinear  problems.  In  contrast,  CEMA  enforces  the  exact  conservation  of  energy  in  the 
nonlinear  regime,  but  at  the  same  time  inherits  high  frequency  dissipation  from  the  underlying  generalized-a 
algorithm.  Consequently,  an  artificial  mechanism  for  transferring  energy  from  the  higher  (artificial)  modes 
to  the  low’er  modes  is  created  that  drives  the  response  to  an  erroneous  solution.  In  contrast,  the  proposed 
ED  scheme  achieves  both  nonlinear  stability  and  high  frequency  dissipation. 

Next,  the  same  problem  was  simulated  with  poo  =  1,  t.e.  with  no  high  frequency  dissipation.  In  this 
case,  two  refinements  in  time  step  size  were  required  to  successfully  complete  the  simulation,  one  at  time 
t  =  0.1665  s  (At  =  5  10“4  s),  the  other  at  time  t  =  0.2142  s  (At  =  2.5  1(T4  s).  Deflections  predicted 
by  the  EP  and  ED  schemes,  shown  in  fig.  5.20,  are  in  good  agreement  during  the  initial  snap-through 
phase,  but  become  increasingly  different  during  the  subsequent  oscillations.  The  force  and  velocity  fields 
are  markedly  different.  The  plate  center  forces  for  both  EP  and  ED  schemes  are  shown  in  fig.  5.21.  The 
forces  predicted  by  the  EP  scheme  present  violent  oscillations  of  amplitude  up  to  an  order  of  magnitude 
larger  than  those  predicted  by  the  ED  scheme.  These  violent  oscillations  hamper  the  convergence  of  the 
Newton  process  at  each  time  step,  leading  to  the  need  for  smaller  time  steps.  The  same  observations  can 
be  made  about  fig.  5.22  which  compares  the  plate  center  vertical  velocity.  Violent  oscillations  are  initiated 
at  snap-through  and  the  strict  preservation  of  energy  implied  by  the  EP  scheme  prevents  any  subsequent 
decay  of  these  vibrations.  Since  vibratory  stresses  are  a  great  importance  to  designers,  it  is  essential  to 
assess  the  ability  of  new  integration  schemes  to  reliably  predict  these  quantities.  It  is  unfortunate  that  many 
scientific  publications  about  geometric  integration  only  present  responses  for  preserved  quantities  such  as 
total  mechanical  energy  or  momentum.  The  above  plots  demonstrate  that  while  EP  scheme  might  perform 
very  well  for  the  prediction  of  total  energy,  momentum,  or  even  displacement  fields,  they  are  unable  to 
reliably  predict  other  important  fields  such  as  velocities  and  internal  stresses.  Consequently,  such  schemes 
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are  of  little  values  in  real  life  applications. 


6.  Conclusions.  In  this  work,  a  robust  algorithm  for  the  dynamic  analysis  of  geometrically  exact  shell 
structures  was  presented.  The  method  is  geometry-based,  i.e.  it  incorporates  knowledge  about  specific 
qualitative  features  of  the  underlying  partial  differential  equations.  However,  departing  from  the  classical 
approaches  based  on  strict  preservation  of  energy,  the  method  presented  here  allows  the  system  to  drift  away 
from  the  level  set  of  constant  energy  in  a  controlled  and  tunable  manner. 

This  feature  achieves  two  goals.  First,  a  bound  is  placed  on  the  total  mechanical  energy  of  the  discrete 
system,  leading  to  the  concept  of  nonlinear  unconditional  stability;  this  stability  criterion  is  stronger  than 
that  obtained  through  the  classical  analysis  of  numerical  schemes.  The  resulting  numerical  procedure  is 
endowed  with  superior  robustness,  an  important  feature  when  dealing  with  complex  engineering  problems. 
Second,  the  monotonic  energy  drift  is  associated  with  numerical  dissipation  of  the  high  frequency  modes. 
This  tunable  dissipation  makes  the  algorithm  stiffly  accurate,  and  avoids  the  build  up  of  energy  in  the  higher 
modes  that  are  an  artifact  of  the  spatial  discretization  process. 

The  proposed  scheme  can  deal  with  general  shell  structures  and  is  not  tied  to  a  specific  spatial  dis¬ 
cretization  of  the  governing  partial  differential  equations.  Kinematic  nonlinearities  are  treated  in  a  rigorous 
manner,  and  material  nonlinearities  can  be  handled  when  the  constitutive  laws  stem  from  the  existence  of  a 
strain  energy  density  function.  The  efficiency  and  robustness  of  the  proposed  approach  were  demonstrated 
with  specific  numerical  examples. 
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Appendix  A.  Rodrigues  Parameters. 

A  common  representation  of  finite  rotations  [21]  is  in  terms  of  Rodrigues  parameters  r  =  2k  tan  0/2, 
where  0  is  the  magnitude  of  the  finite  rotation  and  k  the  components  of  the  unit  vector  about  which  it  takes 
place.  The  following  notation  is  introduced  ro  =  cos2  0/2  =  1  /  (1  +  rT  r/4),  and  the  finite  rotation  tensor 
R  then  writes 


R(r)  -  I  +  r0  r  +  y  rr. 


The  following  decomposition  of  the  rotation  tensor  is  extensively  used  in  this  work 


(A.l) 


(A-2) 


Appendix  B.  Orientation  of  a  Unit  Director. 

Consider  a  unit  vector  £3,  called  a  director ,  that  rotates  to  a  final  orientation  e3.  For  convenience,  this 
director  is  considered  to  be  the  third  unit  vector  of  a  triad  S  defined  by  ix,  i2, 13,  rotating  to  a  triad  S*  with 
orientation  e1}  e2,  e3.  The  relationship  between  these  two  triads  is  ea  =  R  ia,  where  R  is  an  orthogonal 
rotation  tensor.  If  one  solely  focuses  on  the  director,  this  rotation  tensor  is  not  uniquely  defined,  as  any 
rotation  about  the  director  leaves  its  orientation  unchanged.  A  virtual  change  in  the  director  orientation  is 


5e3  =  ef 'Sip, 


where  Sip  is  the  virtual  rotation  vector,  Sip  =  5RRT . 

The  components  of  the  virtual  change  in  director  orientation  measured  in  S*  become 


RTSes  =  RTe%Sip  =  !%RT Sip  =  Stj/ 


■  IT Ao/»*  - 
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(B.l) 


(B.2) 


where  Sip*  are  the  components  of  the  virtual  rotation  vector  in  <S*.  This  relationship  clearly  demonstrates 
that  arbitrary  values  of  <503,  corresponding  to  virtual  rotations  of  the  director  about  its  own  orientation, 
will  not  aifect  virtual  changes  in  the  director  orientation,  and  hence,  setting  Sipl  —  0  is  a  valid  choice.  The 
following  notation  is  adopted 

61P*  =  i±Sal  +i2Sa*2  =  bSa*;  b  =  [il5i2].  (B.3) 

Sa*  is  a  2  x  1,  “two  parameter”  virtual  rotation  vector.  It  follows  that  fop_  =  R  Sip*  =  Rb  So*,  and  hence 

Sea=R7ZbSa*.  (B.4) 

If  Rodrigues  parameters  are  used  to  parameterize  R,  an  equivalent  expression  can  be  obtained  for  finite 
changes  in  director  orientation  with  the  help  of  eq.  (A. 2) 

£af  -Sai  =  Rm’ii  b  s*  =  Qm  s *;  r*  =  b  s* ,  (B.5) 

where  r*  are  the  Rodrigues  parameter  measured  in  S* ,  and  s*  the  corresponding  “two  parameter”  incremental 
rotation  vector. 
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